home *** CD-ROM | disk | FTP | other *** search
/ Mac Magazin/MacEasy 21 / Mac Magazin and MacEasy Magazine CD - Issue 21.iso / Wissenschaft & Technik / yorick12vr1-nofpu folder / include / fitrat.i < prev    next >
Text File  |  1996-02-13  |  5KB  |  160 lines
  1. /*
  2.    FITRAT.I
  3.    Polynomial and rational function interpolators.
  4.  
  5.    $Id$
  6.  */
  7. /*    Copyright (c) 1994.  The Regents of the University of California.
  8.                     All rights reserved.  */
  9.  
  10. func fitpol(y, x, xp, keep=)
  11. /* DOCUMENT yp= fitpol(y, x, xp)
  12.        -or- yp= fitpol(y, x, xp, keep=1)
  13.      is an interpolation routine similar to interp, except that fitpol
  14.      returns the polynomial of degree numberof(X)-1 which passes through
  15.      the given points (X,Y), evaluated at the requested points XP.
  16.  
  17.      The X must either increase or decrease monotonically.
  18.  
  19.      If the KEEP keyword is present and non-zero, the external variable
  20.      yperr will contain a list of error estimates for the returned values
  21.      yp on exit.
  22.  
  23.      The algorithm is taken from Numerical Recipes (Press, et. al.,
  24.      Cambridge University Press, 1988); it is called Neville's algorithm.
  25.      The rational function interpolator fitrat is better for "typical"
  26.      functions.  The Yorick implementaion requires numberof(X)*numberof(XP)
  27.      temporary arrays, so the X and Y arrays should be reasonably small.
  28.  
  29.    SEE ALSO: fitrat, interp
  30.  */
  31. {
  32.   /* find lower and upper indices in x of the interval containing xp */
  33.   nx= numberof(x);
  34.   u= digitize(xp, x);
  35.   l= max(u-1, 1);
  36.   u= min(u, nx);
  37.  
  38.   /* initialize yp to the y(x) closest to xp */
  39.   mask= (abs(x(l)-xp)<=abs(x(u)-xp));
  40.   ns= l*mask + u*(!mask);
  41.   yp= y(ns);
  42.   ns--;
  43.  
  44.   c= d= y= array(double(y), dimsof(xp));
  45.   dx= x - xp(-,);
  46.   if (dimsof(ns)(1)) {
  47.     is= ns;
  48.     is(*)= nx*indgen(0:numberof(xp)-1);
  49.   } else {
  50.     is= 0;
  51.   }
  52.   for (m=1 ; m<nx ; m++) {
  53.     ho= dx(1:nx-m,..);
  54.     hp= dx(1+m:nx,..);
  55.     den= (c(2:nx-m+1,..)-d(1:nx-m,..))/(ho-hp);
  56.     d(1:nx-m,..)= hp*den;
  57.     c(1:nx-m,..)= ho*den;
  58.     mask= (2*ns>=(nx-m));
  59.     y= d(max(ns,1)+is)*mask + c(ns+1+is)*(!mask); /* this is error estimate */
  60.     yp+= y;
  61.     ns-= mask;
  62.   }
  63.  
  64.   if (keep) {
  65.     extern yperr;
  66.     yperr= y;
  67.   }
  68.  
  69.   return yp;
  70. }
  71.  
  72. func fitrat(y, x, xp, keep=)
  73. /* DOCUMENT yp= fitrat(y, x, xp)
  74.        -or- yp= fitrat(y, x, xp, keep=1)
  75.      is an interpolation routine similar to interp, except that fitpol
  76.      returns the diagonal rational function of degree numberof(X)-1 which
  77.      passes through the given points (X,Y), evaluated at the requested
  78.      points XP.  (The numerator and denominator polynomials have equal
  79.      degree, or the denominator has one larger degree.)
  80.  
  81.      The X must either increase or decrease monotonically.  Also, this
  82.      algorithm works by recursion, fitting successively to consecutive
  83.      pairs of points, then consecutive triples, and so forth.
  84.      If there is a pole in any of these fits to subsets, the algorithm
  85.      fails even though the rational function for the final fit is non-
  86.      singular.  In particular, if any of the Y values is zero, the
  87.      algorithm fails, and you should be very wary of lists for which
  88.      Y changes sign.  Note that if numberof(X) is even, the rational
  89.      function is Y-translation invariant, while numberof(X) odd generally
  90.      results in a non-translatable fit (because it decays to y=0).
  91.  
  92.      If the KEEP keyword is present and non-zero, the external variable
  93.      yperr will contain a list of error estimates for the returned values
  94.      yp on exit.
  95.  
  96.      The algorithm is taken from Numerical Recipes (Press, et. al.,
  97.      Cambridge University Press, 1988); it is called the Bulirsch-Stoer
  98.      algorithm.  The Yorick implementaion requires numberof(X)*numberof(XP)
  99.      temporary arrays, so the X and Y arrays should be reasonably small.
  100.  
  101.    SEE ALSO: fitpol, interp
  102.  */
  103. {
  104.   /* find lower and upper indices in x of the interval containing xp */
  105.   nx= numberof(x);
  106.   u= digitize(xp, x);
  107.   l= max(u-1, 1);
  108.   u= min(u, nx);
  109.  
  110.   /* initialize yp to the y(x) closest to xp */
  111.   mask= (abs(x(l)-xp)<=abs(x(u)-xp));
  112.   ns= l*mask + u*(!mask);
  113.   yp= y(ns);
  114.  
  115.   exact= (x(ns)==xp);
  116.   list= where(exact);
  117.   if (numberof(list)) {
  118.     yexact= yp(list);
  119.     yerr= 0.0*yexact;
  120.   }
  121.  
  122.   list= where(!exact);
  123.   if (numberof(list)) {
  124.     yp= yp(list);
  125.     xp= xp(list);
  126.     ns= ns(list);
  127.  
  128.     c= d= y= array(double(y), numberof(xp));
  129.     d+= fitrat_tiny;  /* I question whether this ever rescues a
  130.              failed fit, but I'll leave it here anyway.  */
  131.     dx= x - xp(-,);
  132.     ns--;
  133.     is= nx*indgen(0:numberof(xp)-1);
  134.     for (m=1 ; m<nx ; m++) {
  135.       di= d(1:nx-m,);
  136.       ci= c(2:nx-m+1,);
  137.       t= dx(1:nx-m,)*di/dx(1+m:nx,);
  138.       dd= (ci-di)/(t-ci);
  139.       d(1:nx-m,)= ci*dd;
  140.       c(1:nx-m,)= t*dd;
  141.       mask= (2*ns>=(nx-m));
  142.       y= d(max(ns,1)+is)*mask + c(ns+1+is)*(!mask);
  143.       yp+= y;
  144.       ns-= mask;
  145.     }
  146.  
  147.   } else {
  148.     yp= [];
  149.   }
  150.  
  151.   if (keep) {
  152.     extern yperr;
  153.     yperr= merge(yerr, y, exact);
  154.   }
  155.  
  156.   return merge(yexact, yp, exact);
  157. }
  158.  
  159. fitrat_tiny= 1.e-30;
  160.